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BACKGROUND OF THE INVENTION 



1. Field of the Invention 

[0005] The present invention pertains generally to methods used in the inversion of 
data, more particularly to the inversion of seisniic data, and still more particularly to the 
full wave seismic inversion of seismic data without using a source term. 

2. Description of the relevant art 

[0006] Common seismic industry practice is that velocity structure is estimated by 
analyzing the travel time of a set of seismic signals. In cross hole and surface-to- 
borehole applications, typical approaches involve ray tracing in which the ray may be 
straight or curved depending on the degree of resolution desired, and more recently the 
Fresnel volume approach. The travel time approach is fundamentally a high-frequency 
approximation with a maximum resolution on the order of a wavelength, or a fraction 
(5%) of the well separation in some practical cases. The usefulness of the result obtained 
from the ray tomography may be limited if the objective is to better understand the 
petrophysical and hydrological 'properties' of the soil and rock: an increasingly 
important subject in characterizing petroleum and geothermal reservoirs, with potential 
environmental applications. 

[0007] A better alternative to the travel time approach appears to be full waveform 
inversion. A number of recent studies on this subject suggest that full waveform 
inversion may provide improved resolution of velocity and density structures. Amplitude 
and phase information contained in the waveform are both sensitive to intrinsic energy 
loss of the propagating wave through dispersion, and to the petrophysical properties of 
the material through which the wave propagates. These sensitivities are utilized by full 
waveform analysis to form an attractive tool for investigating hydrological and 
petrophysical properties of a geophysical medium. 

[0008] Full waveform inversion has a major impediment to easy implementation. In 
all field applications, the effective source waveform and the coupling of the medium 
between the source and the receiver are not very well understood. The effective source 
waveform and its coupling to the medium under study are generally taken together as a 
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lumped "source term". The source/medium coupling problem can be resolved to some 
extent by using a good velocity approximation in geophysical applications; however, the 
measured signal cannot generally be accurately calibrated, rendering full waveform 
inversion a technical challenge. 

BRIEF SUMMARY OF THE INVENTION 

[0009] This invention provides for an iterative method of using normalized waveform 
inversion to obtain a model parameter describing one or more physical properties of a 
medium, the method comprising the steps of: a) inputting one or more each of NS source 
and NG measurement spatial locations; b) measuring time domain data at each of the NG 
measurement locations resulting from an input waveform at one of the NS source 
locations propagating through a medium, i) for each of the NS source locations, ii) 
thereby forming a time domain measurement data set D^. (0, 0 = 1 ^G; / = 1 - NS) , 
where d denotes a data tensor; c) Fourier transforming the time domain measurement data 
set D^j.(t) to create a measurement spectral data set D^.(ty) having frequency and 
amplitude information for each of the NG measurement locations; d) normalizing the 
measurement spectral data set D^,.(a^) to create a normalized data wavefield Tj^io)) ; e) 

modeling the medium using an iterated model parameter m describing one or more 
physical properties of the medium, the NS source and the NG measurement spatial 
locations used as respective model input and model response spatial locations contained 
within the model of the medium, and initializing the iterated model parameter m with 
corresponding one or more known bulk properties of the medium being modeled, said 
modeling step comprising: i) creating a measurement model by: (1) applying a delta 
function source collocated with the ^ NS source, (2) modeling the response at the NG 
measurement locations, using a time domain modeling method, to create a synthetic 
medium response at the receiver due to the source, i^)^^ denoting the model 

response, (3) repeating the applying and modeling steps at each of the NS source 
locations and NG measurement locations until the measurement model is full, and ii) 
Fourier transforming the model response (0 ' ob^^J^ ^ frequency domain synthetic 
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response P"(fi^) ; iii) foraiing a normalized modeled wavefield using the frequency 

domain synthetic response T^^ (^) = Pyv (^)[Pu(^)] ^ 5 0 minimizing a weighted 
error, between the normalized data wavefield Tj. (a)) and the normalized modeled 
wavefield T^(fi^) of the response of the medium, to a level below an error bound, i) said 

weighted error met by using the iterated model parameter m , known as the a minimized 
model parameter m . 

[0010] In the method described above the modeled normalized wavefield Tj^(fi^) is 

independent of the input waveform applied to any of the NS source locations, and may be 
used for a variety of applications, so long as they are capable of modeling the physical 
medium under test. Such modeling may utilize finite elements of various types, boundary 
elements, mixtures of both, and/or discrete finite differences. The results of the models at 
the iterative step concluding the minimization results in the minimizing model parameter 
m , which may be directly stored in a computer-readable medium, or immediately 
processed into a display of the minimizing model parameter m as a graphical 
representation of one or more of the physical properties of the medium measured, and 
optionally stored as a data set or data structure in a computer-readable medium. Although 
the thrust of this invention is highlighted herein to the field of seismic analysis, it is not 
restricted to such field. 

[001 1] The methods described above may be further stored as a computer program in 
at least one computer-readable medium, to allow for execution on a computer. 
[0012] The method of normalized waveform inversion uses minimization of 
weighted error (m) that is calculated by a square of the weighing matrix 

multiplied by the difference between the normalized data wavefield T^^ {(o) and the 

normalized modeled wavefield (eo) , or (2> (w ) = || (x; (o)) - Tjj (a)))^ sunmied over 
all frequencies, NS sources, and NG receivers. The weighing matrix is an inverse of 
the standard deviation of the normalized data wavefield TjJ (ty) , or other appropriate 
weighting method indicative of the convergence of the minimization. 
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[0013] The detailed normalization step described above is further comprised of the 
steps of: a) choosing one of the measurement spectral data set D^, (£y) at one of the 

measurement locations, such as Df^ {o)) , i) said chosen measurement spectral data set 

D{^. [q)) known as a reference measurement; ii) dividing the measurement spectral data 

set D^.(6c;) at a particular frequency in the measurement spectral data set D^j(a;) , iii) by 

the reference measurement [co) at the same frequency, b) to create a normalized data 

wavefield Tj.{co). 

[0014] The error bound of the minimizing step is met by first assigning an initial 
weighted error E^, then dividing each subsequent iterated weighted error E^^ by the initial 

E E 
error E^ to create a relative error — . When the relative error — is below a pre- 

estabhshed error bound, preferably below 10*^, more preferably below 10"^ , still more 
preferably below 10'^ , and possibly below 10'^, the iterations of the minimization step 
conclude. Typically an error bound on the order of 10*^ works particularly well. The 
model parameter m that achieved this error bound is then the solution to the 
minimization problem, and the wavefield inversion is complete. 
[0015] In another embodiment of this invention, modeling takes place not in the time 
domain, but rather in the frequency domain as a modal analysis. In this case, the delta 
function input, which has an infinite number of Fourier frequencies, is truncated at the 
highest frequency present in the model data. In this implementation, the frequency 
domain synthetic response PJ(^) is directly obtained in the frequency domain. As with 

the previous methods, the frequency-based model, minimized model parameter 
m solution, and plots of the physical properties of the model parameter m solution may 
be stored on a computer-readable medium for further computations on a computer. 
[0016] In yet another embodiment, an article of manufacture comprises a computer- 
readable medium having the method of full waveform inversion described above placed 
upon it for subsequent execution on a computer. 
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BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE 

DRAWINGS 



[0017] The invention will be more fully understood by reference to the following 
drawings, which are for illustrative purposes only: 

[0018] Fig. 1 A is a 2-D velocity model using a fault model in a background of 3000 
m/s constant velocity 

[0019] Fig. IB is the inversion result using impulse response of the model of Fig. 1 A. 
[0020] Fig. IC is the inversion result using normalized wavefield obtained for the 
model of Fig. lA. 

[0021] Fig. 2 is a root mean square convergence showing the root mean square (rms) 
error and the associated Lagrange multipliers as a function of the number of iterations 
during the full waveform inversion. 

DETAILED DESCRIPTION OF THE PREFERRED 

EMBODIMENT 

Defined terms 

[0022] Computer means any device capable of performing the steps developed in 
this invention to result in an optimal waterflood injection, including but not limited to: a 
microprocessor, a digital state machine, a field progranmiable gate array (FGPA), a 
digital signal processor, a collocated integrated memory system with microprocessor and 
analog or digital output device, a distributed memory system with microprocessor and 
analog or digital output device connected with digital or analog signal protocols. 
[0023] Computer readable media means any source of organized information that 
may be processed by a computer to perform the steps developed in this invention to result 
in an optimal waterflood injection, including but not limited to: a magnetically readable 
storage system; optically readable storage media such as punch cards or printed matter 
readable by direct methods or methods of optical character recognition; other optical 
storage media such as a compact disc (CD), a digital versatile disc (DVD), a rewritable 
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CD and/or DVD; electrically readable media such as programmable read only memories 
(PROMs), electrically erasable programmable read only memories (EEPROMs), field 
progranmiable gate arrays (FGPAs), flash random access memory (flash RAM); and 
remotely transmitted information by electromagnetic or optical methods. 
[0024] Physical property means a derived or observable, measurable characteristic of 
sample, including, but not limited to speed of sound in a specified direction, density, 
permeability, porosity, resistivity, permittivity, Young's modulus, tensor stress-strain 
properties, speed of sound propagation, shear modulus, amount of water present, and 
amount of oil present. 

[0025] Bulk property means a physical property of a medium spatially averaged over 
a volume of interest. 

[0026] Time domain measurement data set means a collection of data measured in 
the time domain such as a seismogram, possibly including spatial information of one or 
more source and/or receiver locations. 

[0027] Wavefield means modeled or measured data, or combination of such data, 
that obey differential equations describing wave motion in time domain or in frequency 
domain. 

[0028] Normalized data wavefield means measured data at any receiver point 
normalized by measured data at a reference point in the frequency domain. 
[0029] Normalized modeled wavefield means modeled (calculated) data at any 
receiver point normalized by modeled data at a reference point in the frequency domain. 

Introduction 

[0030] Source-independent full waveform inversion of seismic data can lead to high- 
resolution imaging of subsurface geophysical parameters that would otherwise have not 
been possible using other conventional methods. Detailed information of distribution of 
geophysical parameters are crucial in understanding the spatial distribution of gas and 
fluid, and their movement when repeated measurements are made over time. Potential 
applications of the proposed method include, but are not limited to exploration of oil, 
methane hydrate, and geothermal resources, characterization and monitoring of various 
petrological reservoirs, CO2 sequestration, and environmental remediation processes. 
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Normalized Wavefield Analysis 



[0031] Without loss of generality, let us assume a field survey involving NS 
transmitter positions and NG receiver positions of arbitrary configurations, but with 
known spatial locations. Additionally, the spatial coordinates of the boundary of the 
medium are also known. In acquiring data, the source and receiver deployment 
configuration is an important part of survey design for successfully achieving the desired 
survey objective. The full wavefomi inversion scheme described herein is general, since 
any arbitrary configuration is acceptable: surface or single borehole reflection, surface-to- 
borehole or borehole-to-surface conventional processing of time-lapse Vertical Seismic 
Profiling (VSP) data, cross-hole, or combinations of these. 

[0032] To develop the full 3-D normalized wavefield inversion method, a full tensor 
measurement is required with three linearly independent component measurements at 
each receiver spatial position for each one of three transmitting events at each source 
spatial position. The data tensor may be in the form of pressure, displacement, velocity 
or acceleration, and may be described in general as 

/>J(0 = /f,(0®P,^(0®5,(0, (j^l-NG; i^l-NS). (1) 

[0033] The superscript d indicates data and each constituent in this equation is a 
(3x3) tensor. Equation (1) states that the data D^^(t) recorded at the receiver position 

is a multiple convolution of the actually transmitted source 5, (0 that includes source 
system function and the radiation pattern caused by source-medium coupling, the impulse 
response Pf^it) of the medium at the receiver due to the source, and the receiver 

system function Rj(t) including the medium-receiver coupling. In the following method 
of analysis, we will remove the Rj(t) by assuming that the receiver (geophone) 

calibration is known and the effect of medium-receiver coupling to data is relatively 
minor compared to that of the source waveform. The source function 5, (0 is generally 
not well understood, and is a subject of considerable ongoing research just by itself, in 
order to improve the quality of data sampling. The impulse response Pf^ (t) is the 
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solution to the governing differential equation of the medium under test, with an impulse 
source in time at the i"* position. If we Fourier transform Equation (1), 
FT{{DJSJ*){t)} (D,S,P)(<y) , with the Rjit) term neglected, we get 

D;,(a;) = Pj(^y)S,(fi;). (2) 

[0034] The transformed frequency domain data contains all of the information the 
time series originally contained. In the time-series data, longitudinal and transverse 
waves of primary arrivals are followed by events of converted modes, multiple 
reflections, and all other scattering waves caused by heterogeneities in the sampling 
medium. The data at the receiver position may then be described with Equation (2) as 

I>?.M = K.<,dy(a;) = 



di.. 



\cji 



(3) 



Multiple subscripts HkjV are used to describe data elements. The last two subscripts 'jV 
only indicate receiver (/ = 1, 2,..., NG) and source {i = 1, 2,..., iV5) positions, and the first 
two subscripts Hk' are used to describe Cartesian field components / = (1, 2, 3). Elements 
of the tensor source function S^{w) in equation (2) are unknown. The tensor source 

function is composed of three source vectors [s^pS^pS^J(ty) and each source vector may 

in turn be composed of three unknown Cartesian components; s^, = {sxa^s^^^^s^^^ , 
So, the source tensor at the i-ih position may be fully described by 



S,(^) = {s.,>s,,,s,J(fi;) = 



^\ai ^\bi ^ici 
^2ai ^2bi ^2ci 



(4) 



[0035] The tensor impulse response of the medium that relates the diagonal impulse 
source at the i* transmitter position to the measurements at the f* receiver position may 
be written as 
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Puji Pnji P\3ji 



}(^)= P2.;. Ptlji P^Ji (^)- 



(5) 



Pyiji Pyiji P33ji 



[0036] To define the normalized wavefield (alternatively known as the transfer 
function of the medium under test), we first select the reference point (which may be any 
of the measurement locations), say 7 = 1 for convenience. The tensor normalized data 
wavefield (or data transfer function) is now defined by the data aty = 1 - NG, normalized 
by that of the reference point: 



where the convention of the subscripts for the transfer function is the same as that used 



transfer function is defined as the normalized impulse response of the medium; hence the 
uniqueness of the transfer function. The necessary condition for the source term to divide 
out is that the determinant of the source matrix is non-zero. In other words, the source 
events need to be linearly independent. This condition can be met, in principle, even if 
source events consist of explosions, so long as their Cartesian constituents are linearly 
independent. 

[0037] In this procedure we assume that for each source, / = 1-- NS, NG 
measurements are made simultaneously. If for some reason all NG data cannot be taken 
simultaneously, a simultaneous data set may be simulated if neighboring subsets share at 
least one overlapping receiver position. Since we do not know the medium, the impulse 
response itself is not known either at this point. 



-Pj(M^)S,(M^)[S,(iv)r[P.f(M^)] 



1-1 




(6) 



for impulse response. In Equation (6) the source term S, {co) divides out to unity, so the 
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Full Waveform Inversion Using Normalized Wavefield Analysis 



[0038] To obtain the numerical solution for the media impulse response for scalar 
and elastic wave equations one needs to spatially discretize the constitutive parameters, 
and apply finite difference, finite element, or integral equation (boundary element) 
techniques to solve the governing differential equation of the discretized system. Using 
the finite difference or finite element method on a mass/spring/damper model of the 
medium under test, the assembled system of equations, including the damping, takes a 
general form: 



where the field vector p[t) is the discretized wavefield, M is the mass matrix, K is the 

stiffness matrix, and C is the damping matrix. If there are a total of N unknowns in the 
discretization, then all M, C, and K matrices are NxN square, the field vector/? is Nxl, 
and the load vector g is also ^xl whose entries are all zero except for T values at the 
source locations. Boundary conditions are implicitly included by the known spatial 
coordinates, and any level of continuity of the derivatives of the boundary using well- 
known finite element, boundary element, or finite difference methods, depending on the 
model solution method utilized. 

[0039] The reduced system of equations may be solved in the time domain, typically 
by using coupled first-order differential equations on a staggered grid, or in the frequency 
domain after Fourier transforming Equation (7), FT{p(t)} ->p{o)) , into 



[0040] Next, we will show that the normalized wavefield defined by Equation (6) is 
all that is required for full waveform inversion. In the inversion the objective error 
functional typically consists of a measure of the error between the measured data and the 
resulting model fit. In this case, the normalized wavefield can be used to evaluate the 
error term. For a given model one can generate synthetic data using Equation (8), and 



Mp(t)^Cp{t)^Kp{t)^gS{t) 



(7) 



'(o^Mp {co) + iwCp (o)) + Kp (^) = g . 



(8) 
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then obtain the model normalized wavefield similar to the data normalized wavefield 
described by Equation (6). 

[0041] Formally, the synthetic impulse response at the receiver due to the 
source may be obtained and designated as FHico) , 

where the superscript m indicates model. Accordingly, the model normalized wavefield 
(or transfer function) is obtained from the numerical solution for the given velocity model 

{o)). (10) 

[0042] The inversion procedure begins with a measure of the error between the 
normalized wavefield (or transfer function) model and the data T'' , with the 
subscripts to the normalized wavefield dropped 

0{m) = l^N,(T^-r)f. (11) 
The misfit error between data and model transfer functions at the reference point is 
always zero (r^^ii ^^iku = O) . In setting up the error functional of Equation (1 1), both real 

and imaginary parts are separated, so the actual number of data parts used for the 
inversion is NEQ = 2xNFREQxNTx{NG -l)><9 , and the computation is done in real, or 
extended real arithmetic. The variable NFREQ is the total number of frequencies. The 
matrix is an NEQxNEQ weighing matrix that is the inverse of the standard 

deviation of the normahzed data wavefield T'' . 

[0043] For inversion we consider Gauss-Newton method by first expanding the 
objective functional, Equation (11), into a Taylor series 



r% T% n 

Pnji Pi2ji P\$ji 

Pziji Pllji Pl^ji 

P^iji Pyiji Pz^ji 



(9) 



*21/i 

A.fn 

hxji hiji hzji 



13;/ 
m 

23 ji 
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ifi{m'^Sm) = ip{m)'^yJSm'\-0.5Sm''HJm+o{{Smf] . (12) 
Here, is a perturbation to the minimizing model parameter m , Ym ^ 



column matrix consisting of elements 



dm. 



p = 1 - Af > with M being the total number 



of parameters to be determined, and is compactly written as = 2y\V/\V^ (T"* -T*' ) , 
and is an MxM square matrix (a weighted and constrained Hessian) consisting of 



elements 



written as = 2J^ W/ W^J + o\^^ 



[0044] The last term above, O 



f BJ^ 

— , is literally interpreted as the change in the 
\omJ 



partial derivative of data (transfer function in this case) due to a change in the velocity. 
This term is small if either the residuals are small, or the forward differential equation is 
quasi linear. Furthermore, it is usually difficult to compute, so is generally dropped. For 
each frequency and source the sensitivity function (Jacobian) J is a {2 x (NG-l) x 9}x M 
rectangular matrix: {2 x (NG-l) x 9} because real and imaginary parts have been 
separated for each of 9 tensor elements. For example, for the source and a fixed 
frequency, the entries to the Jacobian corresponding to the receiver and the 
parameter may be obtained as 



with 



real 
imaginary; 



dm. 



*'lkji » * , A. — 1, ^, J , 



(13) 



d'", IPw 



I_ tit 
Pll 



dm. 



IPU 



dm 



1 J 



(14) 



1= 1,2, ...,NS, 7 = 2,3, ...,iVG, 9= 1,2, ...,M. 



[0045] Here, | | is the determinant of the tensor impulse response at the reference 
receiver point and { , /, = 1, 2, 3 } are the result of multiplication of p J, {o)) and the 
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adjoint of p'Hia)) (see Appendix A). The partial derivatives of the normalized wavefield 

can be explicitly written as a function of the partial derivatives of the impulse responses 
(see Appendix B). Partial derivatives of the impulse responses may be obtained using the 
forward modeling results using equation (8). The full waveform inversion of seismic 
data based on equations (11) through (14) will not require the knowledge of the actual 
source waveform, which is an advantageous feature of this method. 
[0046] The error functional that will be minimized consists of the error term, 
Equation (12), and a constraint that will have a smoothing effect on the variation of the 
model in the updating process. Specifically, it may be written as 



where A is the Lagrange multiplier that controls relative importance of data errors and 
the behavior of the parameter variation, and W„ , which is an MxM controlling matrix. 

When the W„ matrix is diagonal it has an effect of keeping the parameter from changing 

from the current one. On the other hand, if the matrix represents a gradient operator its 
effect is to spatially smooth out the changes. Minimization of the functional of Equation 
(15) with respect to perturbation in model parameter results in a system of normal 
equations 



from which the model parameter at the (^+l)-th iteration is updated to m**" = w * + Sm . 

[0047] The iteration process stops when the change in the model parameter Sm is 
below a preset tolerance, typically given in terms of root mean square (rms). 

Construction of the Source Function 

[0048] Once the inversion is successful, the functional description of the actual 
source waveform can be established by deconvolution. The deconvolution can be 
achieved by dividing the data with the impulse response of the synthetic model that has 




(15) 



(rw/w,j+ 



(16) 
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been obtained through the inversion process. From Equation (2), the deconvolution 
tensor may be estimated to be 



C,,(^) = Di(ty)[p;]"'(^) = 



"Ibji 



^\cji 



^laji ^2bji ^2cji 
^3aji ^3bji ^3cji 



(17) 



[0049] Notice that the elements of the impulse response of the model Py7 {(o) are 
independent of the source events (a, b, c), but the elements of the deconvolution tensor 
Cy,. (ty) are related to the source events. If the data and the inversion process were exact, 

the deconvolution function Cy, [o)) is actually the frequency spectrum of the source, and 

it must therefore be identical for all receiver position. Since neither the data nor the 
inversion process are exact, the source function may be obtained by taking a statistical 
mean of the deconvolution function at different receiver positions. For example the first 
element of the source function in Equation (4) may be estimated using the root mean 
square approach 



(18) 



[0050] Since the function C^, {(o) is itself complex, the averaging process must be done 
in amplitude and phase, or real and imaginary parts, separately. Inverse Fourier 
transforming, FT~^ {S, {o))] — > 5, (/) , one obtains the effective source function in time: 



S,(O = K,^w>*«}(0 = 



^2ai ^2bi ^2ci 



^3*1 hci J 



(19) 
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[0051] Note that each one of three events {a,b, c) , for example 

^ai {^)^{Ki>hai^haiY (0 » ^^^^ produce a complete description of the effective source at 
the transmitter location / as a function of time. The function (r) describes the source at 
the transmitter position that would most likely have produced data/J^. (/) , y = 1 ~ NG , 

shown in Equation (1). The reconstructed source function is the effective source at the 
time of the event, and will also contain the spatial and temporal radiation patterns. 

Full Waveform Inversion of 2-D Acoustic Velocitv 

[0052] The inversion method described here has been tested using a simple 2-D 
acoustic model. Let us consider the impulse response governed by a 2-D acoustic wave 
equation in the frequency domain 

V^p(x,x,,Q))+^TTPix,x,,Q)) + S(x-x,) = 0, (20) 

V \X) 

where the impulse response p is the scalar pressure wavefield, v the velocity, and (jc,jcJ 
are the field and source positions in 2-D. The source is an impulse point source expressed 
as a 2-D spatial delta function S[x-x^) located at Xs, and is also a delta function at r = 0 

in the time-domain formulation, A finite-element modeling (FEM) scheme is used to 
generate the synthetic impulse response. The model parameter is the acoustic velocity in 
each of the rectangular elements used for the FEM solution. Following the procedure 
described in the previous section, the scalar synthetic normalized wavefield is obtained 
from the numerical solution for the given velocity model 

^;v(^) = '^^> (; = 1-A^G; i^l-NS). (21) 
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The inversion procedure starts with an objective functional, Equation (15), reduced to 
handle scalar problem. The number of equations is NEQ = 2 x NFREQ xNTx (NG - 1) , 
with all computations done in real arithmetic. The related sensitivity functions are 



[0053] Figures 1 A shows a model 100 used for the test. The model is a broken 
dipping fault in a background of material 105 having a 3000 m/s constant velocity. The 
fault consists of a 6 m thick low velocity (2500 m/s) layer 130 overlain by another 6 m 
thick high velocity (3500 m/s) layer 140 starting an x position of 30-60 meters, and 
discontinuously dipping respectively toward the x = - 55 meters position with low 
velocity section 135 and high velocity overlay 145. A cross-hole configuration is used 
for the exercise, with the borehole at a: = - 45 m for the transmitter (Tx) borehole 1 10 and 
the other at a: = 45 m for the receiver (Rx) borehole 120. A total of 21 line sources 
starting at elevation z = 0 are used with an equal vertical separation of 9 m, and the same 
number of receivers and separations for the receivers. Each of the source pressure 
wavefields computed at 21 receiver positions have been normalized by the first pressure 
wavefield, resulting in 21 transfer functions. The number of frequencies used is 10; 
starting from 10 Hz to 100 Hz, linearly separated by 10 Hz. The shortest wavelength in 
the background is 30 m, so the maximum resolution is expected to be on the order of 
7-8 m if we take into account the wavelength as a measure of resolution. 
[0054] Using the numerical solution of the model as synthetic data, the inversion was 
initiated with an initial model of 2850 m/s uniform whole space. A grid of uniform cell 
size, 3 m by 3 m, was used throughout. The inversion domain was 120 m by 180 m, 
containing a total of 2400 velocity parameters. 

[0055] The size of the matrix from Equation (15) is relatively modest for the test 
model, so it was solved by using QR decomposition with successive Householder 
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transformations. The Lagrange multiplier X is automatically selected in the inversion 
process. It starts with executing a given number, say nl, of inversions using nl different 
multipliers that are spaced appropriately. The same Jacobian is used at this step. As a 
result nl updated parameter sets are produced, followed by nl forward model calculations 
resulting in nl error terms. Among these, we choose the model and parameter X giving 
the lowest error term result. 

[0056] To demonstrate the validity of the inversion method disclosed here, we first 
carried out conventional inversion using impulse response, with the result is shown in 
Figure IB, which uses the same numbering scheme of Figure lA. Additional features of a 
slightly slower velocity 150 and a slightly higher velocity 160 are pointed out as artifacts 
of the limited data set reconstruction. Using standard conventions, we assume that the 
source waveform is known, so that the data can be reduced to an impulse response, with 
an inversion carried out by using the impulse response. 

[0057] It should be pointed out that this numerical simulation has an ideal impulse 
response, with ideal coupling of the impulse to the medium. This situation is far from 
conmion in the field, however, does provide a comparison with the transfer function 
methods of this invention, which does not depend on source-to-medium coupling. 
[0058] We additionally obtained the velocity structure by using the transfer function 
approach, with the result shown in Figure IC, which also uses the same numbering 
scheme of Figure lA. Additional features of Figure IB showing a slightly slower 
velocity 150 is no longer readily apparent in the transfer function approach of Figure IC. 
However, the slightly higher velocity 160 of Figure IB is still somewhat apparent, but 
reduced in Figure IC. Again, such artifacts are due to the limited data set reconstruction. 
In this exercise we used nl =3 in each iteration to select parameter update and Lagrange 
multiplier. After 5 iterations, two of the results appear almost identical. Slight 
differences may have been due to the fact that the transfer function approach has one less 
data sample than the impulse response approach because one data sample was used to 
normalize the others. The behaviors of the rms errors are also similar as discussed below. 
[0059] Figure 2 compares in chart 200 the rms errors of both the impulse response 
210 to the transfer function method 220. Note that the Lagrange multiplier has also 
changed as with iteration. Coordinate points 211, 212, 213, 214, and 215 are the first five 
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iterations of the impulse response 210. Similarly, coordinate points 221, 222, 223, 224 
and 225 are the first five iterations of the transfer method 220 described herein. It is seen 
that the transfer method generally yields a better result, and is quite similar in response to 
the ideal impulse response 210 method. The quality of inversion could be improved by 
using higher frequency data and by using a higher density deployment of transmitters and 
receivers, but the method has been successfully demonstrated. 

Conclusions 

[0060] An innovative, rigorous full waveform inversion method has been disclosed, 
and the validity of the method successfully demonstrated using a simple 2-D synthetic 
model. An important advantage of the invented method is that full waveform inversion 
of seismic data may be accomplished without the source information by using the 
properties of the normalized wavefield describing the system (typically a geological 
sample) under test. Under the proper combination of sources and receivers, the 
normalized wavefield is shown to be uniquely determined in terms of the normalized 
impulse response. The methodology is readily extended to include general 3-D problems. 
As an important byproduct, it has also been shown that the source function can be 
reconstructed once the full waveform inversion is completed. The reconstructed source 
function describes the effective source, not the source system output, at the time of event, 
including the spatial and temporal radiation patterns. 

[0061] All publications, patents, and patent applications mentioned in this 
specification are herein incorporated by reference to the same extent as if each individual 
publication or patent application were each specifically and individually indicated to be 
incorporated by reference. 

[0062] The description given here, and best modes of operation of the invention, are 
not intended to limit the scope of the invention. Many modifications, alternative 
constructions, and equivalents may be employed without departing from the scope and 
spirit of the invention. 
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APPENDIX A - The Normalized Wavefield 



[0063] The relationship between the normalized wavefield and the impulse 

response is derived. We first rewrite the normalized wavefield defined in equation (6), or 
equation (10). After dropping the superscript dy the source position index /, and the 
angular frequency O) , we write 











T,=P>[P.]" = 

















(A-1) 



along with the impulse response, equation (5), as 
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[0064] The matrix operation in equation (A-1) may be expanded to 



T;=P;[P.r' = 
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where the determinant of the reference impulse response is given by 
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= P1I1P221P331 + P121P23IP311 + AJIP321P211 

" P3nP22\Pm ~ P32iP23\P\n " PmPnxPiu 
and elements of the reduced tensor normalized wavefield are 
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APPENDIX B - The Partial Derivatives of the Normalized Wavefield 



[0065] In Appendix A, each element of the normalized wavefield was shown to 

be a rational function of combination of elementary impulse responses. Let us rewrite 
equation (14) for the partial derivatives of the elementary normalized wavefield with 
respect to the parameter ntq. After dropping the superscript m and the source position 
index /, we get 
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where the determinant is given by equation (A-4). To illustrate the process of obtaining 
partial derivatives of normalized wavefields, we will choose an element, say , which 

is given in equation (A-5). Partial derivatives of |Pj and with respect to the 

parameter are: 

dm, dm, dm, dm. 
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[0066] Substituting (B-2) and (B-3) into (B-1), we get the partial derivative of r,,^ 

with respect to the parameter m,. So, partial derivatives of the normalized wavefield may 
be evaluated in terms of a combination of partial derivatives of the impulse responses: 



dm^ 
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where NG is the number of receiver measurement positions. 
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